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Crystalline membranes at finite temperatures have an anomalous behavior of the bending rigidity 
that makes them more rigid in the long wavelength limit. This issue is particularly relevant for 
applications of graphene in nano- and micro-electromechanical systems. We calculate numerically 
the height-height correlation function G(q) of crystalline two-dimensional membranes, determining 
the renormalized bending rigidity, in the range of wavevectors q from 10 -7 A^ 1 till 10 A" 1 in the self- 
consistent screening approximation (SCSA). For parameters appropriate to graphene, the calculated 
correlation function agrees reasonably with the results of atomistic Monte Carlo simulations for this 
material within the range of q from 10" 2 A -1 till 1 A -1 . In the limit q -> our data for the 
exponent r\ of the renormalized bending rigidity Kjj(g) oc q~ n is compatible with the previously 
known analytical results for the SCSA r\ ~ 0.82. However, this limit appears to be reached only 
for q < 10~ 5 A -1 whereas at intermediate q the behavior of G(q) cannot be described by a single 
exponent. 

PACS numbers: 81.05.ue, 68.60.Dv, 63.20.Ry, 46.70.Hg 



I. INTRODUCTION 

A very active field in statistical mechanics and con- 
densed matter physics is the study of interfaces and mem- 
branes. Physical membranes are two-dimensional sur- 
faces embedded in three-dimensional space. In these sys- 
tems, the interplay between the two-dimensional geome- 
try and thermal fluctuations is at the origin of a number 
of unexpected behaviors, going from flat to glassy and 
tubular phasesP The stability of a flat 2D phase seems to 
be in contradiction with the Mermin- Wagner theorem P 
which states the impossibility of long range order in two 
dimensions due to thermal fluctuations. This apparent 
contradiction became subject of great interest after the 
discovery of graphene, a single atom thick layer of car- 
bon atomspEl which can be considered as the prototype 
of crystalline membranes. The stability of this material 
against crumpling, demonstrated even for free-standing 
samplesPsl wasproven to be related to the presence of 
intrinsic ripples. 10 Ripples and the mechanical properties 
of graphene have been subject of much recent theoretical 
workPJflZI 

The first attempt to study the anomalous elastic- 
ity in polymerized membranes was done by Nelson and 
Peliti using a simple one-loop self-consistent theory, with- 
out including any renormalization of the in-plane Lame 
constants.^ They found an anomalous bending energy 
of the flat phase that for small wave vectors q deviates 
from its constant value and acquires a power-law behav- 
ior for the effective bending rigidity k_r(q) ~ q~^ with 
rj = 1. The existence of anomalous elasticity was con- 
firmed by an e = 4 — D expansion, where D is the mem- 
brane dimension PS A step beyond was done by Le Dous- 
sal and Radzihovsky^ who generalized to polymerized 
membranes the self-consistent screening approximation 
(SCSA) introduced by Bray^S to estimate the critical ex- 



ponents of the O(n) model in the large-n limit. This ap- 
proximation is exact when the co-dimension d c = d — D 
is going to infinity (d being the dimension of the embed- 
ding space) . In Ref . [20] an approximate solution of the 
SCSA in the long wavelength limit was found, giving an 
exponent rj &s 0.821 for a 2D membrane in a 31? space. 

Motivated by the relevance for graphene, several works 
have recently appeared studying the bending rigidity 
properties of 2D crystalline membranes. Mariani and 
von Oppen studied the one-loop correction to the bend- 
ing rigidity due to the effective interaction between flex- 
ural phononsP^l More sophisticated methods as non- 
perturbative renormalization group (NPRG) have been 
used by Kownack and Mouhanna, who found an expo- 
nent of rj s» 0.85p2l in good agreement with the SCSA 
results p3 and by Braghin and Hasselmann, who extended 
the analysis of Ref. [23] to finite momenta. Further- 
more, the validity of SCSA has been recently checked by 
GazitJ^ who has applied the approximation to second 
order expansion in l/d c and found no significant devia- 
tion from the first order expansion. As a result, vertex 
corrections can be neglected during the calculation and 
SCSA seems to be applicable to crystalline membranes. 

In this paper, we solve numerically the SCSA equa- 
tions for the height-height correlation function G(q) and 
calculate it in a wide range of wavevectors q. In the 
long wavelengths limit q — > 0, our results for the ex- 
ponent r) agree with the analytical solution of Le Dous- 
sal and RadzihovskjE^l but at larger q the full solution 
has a more complex form that cannot be described by 
a single exponent. Furthermore, we identify the length- 
scale separating the harmonic behavior in the short wave- 
length limit, from the region where anharmonic coupling 
start to play an important role and the correlation func- 
tion G(q) is renormalized. We also compare the results 
of the numerical solution to Monte Carlo simulations of 
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graphene based on the LCBOPII bond order potential.^ 
The two approximations reasonably agree, justifying the 
use of SCSA in the calculation of physical properties of 
graphene. 



II. METHOD 

In this section wc briefly review the SCSA for 
membranes! 20 * 24 ^ In the Monge representation, displace- 
ments of a D-dimensional membrane embedded in 
a d-dimensional space, are parametrized using a D- 
component phonon field u, and the out-of-plane height 
fluctuations by a d c = d — D dimensional field h. There- 
fore, if r describes the position of a particle on the 
undistorted (flat) membrane, its configuration after the 
displacement due to perturbations will be given by the d- 
dimensional vector r = (ro +u, h). Assuming an asymp- 
totically flat geometry with small out-of-plane fluctua- 
tions, such that u and h are functions of ro, the free 
energy takes the form^ 



F[u,h] = 



f 



d u v 



[ K (V 2 h) 2 



(1) 



where the strain tensor u a 
dients of u and h, reads 



1 



to the lowest order in gra- 



u ap ~ ^{daUfs + d p u a + d a h ■ dph), 



(2) 



with a, j3 = 1, D. In Eq. Q, k, A and /i are the bend- 
ing rigidity, the first Lame constant, and the shear mod- 
ulus respectively.^ In the harmonic approximation, the 
last term of Eq. (p| is neglected, leading to a decoupling 
of the bending (h) and stretching (u) modes. Eq. ([!]) 
provides a correct description of elastic free energy and 
height fluctuations of a membrane as long as the equilib- 
rium phase is truly a flat phase. Once the phonons have 
been integrated out, the effective free energy can be ex- 
pressed in terms of the Fourier components of the height 
fields 



FeffM = 



d D l 

(27T) 



nq 
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d D k' 



i? (Z3) (k,k',q)(h k • h q _ k )(h.k' • h q k , 



(3) 



where the effective four-point-coupling fourth-order ten- 
sor R( D \k, k', q) reads 

J R^(k,k',q)=2 A1 [kP T (q)k'] 2 



+ 



2^A 
2(i + A 



[kP T (q)k][k'P T (q)k'], (4) 



and -Pj^(q) = (S a p — q a qp/q 2 ) is the transverse projec- 
tion operator. Notice that the interaction is completely 
separable for physical membranes (D — 2 and d — 3), al- 
lowing us to writePi?( 2 )(k,k',q) = 26 [q x k] 2 [q x k'] 2 , 
where q = q/q and b = 2(i(fi + A)/(2// + A). 
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FIG. 1. (Color online) Integrand of Eq. Q for p = 2.5 A -1 
and <j> = 0. The UV cutoff in this figure has been taken, for 
illustrative reasons, to q ma x = 5 A -1 . 



Our aim is to calculate the correlation function 



(h a (-q)h f3 (q)) = S a pG(q), 



(5) 



nq +E(q), where S(q) is the self-energy 
4 is the correlation function in the har- 



with G _1 (q) 
and Gq 1 (q) = nq 
monic approximation. In the SCSA theory, the renormal- 
ized elasticity is determined through a 1 /<i c -expansion for 
the 2-point and 4-point correlation functions of h, that 
turns them into a closed self-consistent set of coupled in- 
tegral equations for the self-energy S(q). For physical 
membranes, the set of equations can be written asP3 



G- 1 (q) = Go 1 (q) + S(q) 



(6) 



S(q) 
6(P) 
Hp) 



d 2 p 

J2nf 
bo 



1 + 36 I(p) 
1 f d 2 q 



Kp)[q^ T (p)q] 2 G(q- P ) (7) 

(8) 



(27T) 



g1p-qrG(q)G(p-q) (9) 



In Eq. ([8| the constants k, A and /i appearing in bo are 
divided by ksT, where T is the temperature and fee the 
Boltzman constant. Eqs. ([6|-([9]) admit an analytic solu- 
tion in the long wavelength limit, under the assumptions 
G^ 1 (q) ps T,(q) ps Z/q 4 ~ n , with Z a non-universal ampli- 
tude, and b(k) ps l/3J(fc). The solution of such simpli- 
fied system gives for the critical exponent r\ — 0.82lP3 
However, a full knowledge of the correlation function is 
lacking in this approach. 



III. RESULTS AND DISCUSSION 

In the following, we solve numerically the set of equa- 
tions Eq. Q-Q. The self-consistent cycle starts with 
the harmonic approximation G(q) = Go(q). From this, 
we compute Eq. (JtJ) — <j9j) and the obtained self-energy 
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FIG. 2. (Color online). Evolution of the calculated Gi(q) 
for each iteration i = 1, 51 (red lines). Go(q) = n/q 4 is 
denoted by the dotted blue line. Inset: Go (go) — Gi(qo) as a 
function of the iteration i for qo = 10 -7 A -1 , which shows 
how the solution converges after a few iterations. 



S(q) is used to dress the new correlation function G(q), 
which in turn allows us to start a new iteration. Tak- 
ing into account that G(q), £(q), 6(q) and /(q) depend 
only on the modules of the vector variables, it is nat- 
ural to integrate Eqs. ([6|-([9| in polar coordinates with 
the replacements p — > (p, <p) and q — >• (q,ip). Moreover, 
further in this paper we will make no difference between 
G(q) and G(q). Thus, Eqs. ([6])-(|9]) can be written as 
follows: 



£(g,V0 



G- 1 (q,^) = Go 1 ( g » + S(9,^) 
1 



(10) 



27T (-IJmax 

d0 / dpb(p, (j))pq A sm\ip — (f>) 
o Jo 

xG(\/g 2 +p 2 - 2qpcos(ifj -4>),ijj- <p) (11) 

b (P^)= ^J° T( (12) 
1 + 3M (p, <p) 



—2 / #/ dgg 3 (<z 2 +p 2 -2 M cos(0-^)) 



xG(q,^)G(^q 2 +p 2 -2pq 



cos 



V0,<£-V0 (13) 



In the numerical implementation we have used a (hard) 
ultraviolet (UV) cutoff <7 max in the radial integrals. Due 
to finite size effects, it is natural to consider an UV 
cutoff (which is of the order of the inverse lattice con- 
stant in crystalline membranes) and we have checked 
that the results are independent on this cutoff. We have 
checked that, in the relevant range, the same results 



FIG. 3. (Color online). Comparison of the unrenormalized 
correlation function in the harmonic approximation Go(q) 
(dotted blue line) to the solution of the SCSA equations G(q) 
(red line) and the long- wavelength limit solution G apP (q), 
using the approximations of Ref. (g ra y line). The 

black dashed line is a fitting to the approximate solution 
G(q) w Z/q 4 ~ r ' choosing the parameters r\ = 0.821 and 
Z = 1.2. The vertical dot-dashed line indicates the wavevec- 
tor q* ~ 0.18 A -1 obtained from the Ginzburg criterion Eq. 



are obtained by multiplying G(q) by a cutoff function 
A(q) ~ e - K (i/ A ) j w here A ~ QWx/S an d we have used 
q max = 100 A" 1 . 

The next difficulty is the divergence of the correlation 
function G(q) in the infrared (IR) limit when q — > 0. 
As an example, Fig. [I] shows the integrand of Eq. |9j 
where the two divergences, for q — and q — p, can 
be seen. To avoid such IR divergence, we replace the 
function Go(q) — n/q 4 by Go(q) = n/(q + e) 4 , where e 
is a small parameter (e = 1CP 46 A" 1 in our numerical 
calculations) . 

Because of the power law behavior of the correlation 
function, it is extremely convenient to use a logarithmic 
grid for numerical evaluations. Therefore we discretize 
the momentum axis into points yj = ae d ^ i_1 \ where i is 
the index of the point in the grid of q, a is the minimum 
value considered for the representation (a = 1CP 7 A -1 
in our calculations) and I — log(q max /a) / (N — 1), where 
<7max is the UV cutoff and N is the number of points in the 
grid ofgP3lnFig.|] we show the renormalized correlation 
function G(q) after each of the first 51 iterations. In 
general, convergence is very fast and achieved after about 
20 iterations. 

Our results are summarized in Fig. [3] There we 
compare the bare (unrenormalized) correlation function 
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Go(q) = 1/nq 4 (dotted blue line) to the solution of the 
SCSA, G(q) (red line). The important result is the value 
of the wavevector, q c w 0.1 A -1 , where G{q) changes be- 
havior from harmonic, where G(q) oc 1/g 4 (for q > q c ), 
to non-harmonic, with G(q) cx l/q 4 ~ n , for q < q c . The 
Ginzburg criterion^ gives an approximate value of the 
wavevector q* , and thus the spatial scale, L* = 2n/q* , at 
which anharmonic effects become dominant 



q = 



3TK 

87TK 2 ' 



(14) 



where K is the 2D bulk modulus. For graphene, K = 
12.4 eV • A -2 and k = 1.1 eV at room temperature (T = 
300 K)P! leading to q* « 0.18 A" 1 . This wavevector is 
represented by the vertical dotted-dashed line in Fig. [3j 
and it is in good agreement with the SCSA results. 

Furthermore, we have numerically solved the SCSA 
set of equations Eqs. ^ (|9| in the long wavelength ap- 
proximation used by Le Doussal and RadzihovskyPil By 
taking G^ 1 (q) « £(g) and b(p) « l/3/(p), we obtain the 
approximate solution shown by the green line in Fig. |3j 
which is only valid in the long wavelength limit. No- 
tice that both, the exact and the approximate solutions 
coincide for small wavevectors (i.e. in the limit q — > 0). 

Finally, we have fitted this approximate solution to 
G(q) « Z/q 4 -^, with 77 = 0.821 and Z = 1.2, as shown by 
the dashed black line (q is expressed in A -1 ). The three 
results (exact numerical solution of the SCSA, approxi- 
mate numerical solution and analytic approximation) co- 
incides in the long wavelength limit, and corroborate the 
value given in Ref. [2111 for the critical exponent, rj — 0.821. 
We mention here that the above solution is robust as far 
as we start the first iteration from the harmonic approx- 
imation [Go (3) ~ <7 4 ] or from any correlation function 
that diverges faster than q~ 4+r io with 770 ~ 0.85. 

We also compare the solution of the SCSA system of 
equations with the correlation function G(q) of graphene 
extracted from the Monte Carlo simulations presented in 
Ref. 03] For more details about the Monte Carlo cal- 
culation of the correlation function G(q), see Ref. [28] 
In Ref. Q3I the results for the correlation function found 
for two different model potentials were described by a 
power law with exponent 77 = 0.85. The Monte Carlo re- 
sults are shown in Fig.|3]together with the solution of the 
SCSA system of equations and the unrenormalized corre- 
lation function Go(q) = l/nq 4 . In Fig. [4] we can see that 
G(q) obtained from the SCSA equations agrees rather 
well with the Monte Carlo data in the range of q accessi- 
ble in atomistic calculations. An even better agreement 
with Mote Carlo data was found in Ref. [T6] where the 
height-height correlation function was computed using a 
more accurate approximation as the NPRG. However, 
notice that we do not use here any additional adjustable 
parameter when comparing to Monte Carlo data. There- 
fore, this justify the use of SCSA in the intermediate 
range of momenta. 

Furthermore, we compare the results to the approxi- 
mate correlation function G a (q), obtained from the effec- 
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FIG. 4. (Color online), (a) Comparison of the unrenor- 
malized correlation function in the harmonic approximation 
Go{q) (dotted blue line) to the solution of the SCSA equations 
G(q) (red line) and the Monte Carlo data (black dot-dashed 
line). The dashed green line correspond to the approximation 
given by Eq. ( 15 1. In the inset we show the deviation of the 



approximation G a (q) from the SCSA solution G(q). (b) Zoom 
of Fig. 4(a) focusing on the comparison of G(q) from SCSA 
to the Monte Carlo data. 



tive Dyson equation^ 

G- 1 (q) = Gv 1 (q) + X(q), 



(15) 



where Go(q) is the correlation function in the harmonic 
approximation 



Go(q) 



N 



nS q 4 



(16) 
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N being the number of atoms of the sample and So = 
L x L y /N the area per atom, and the self-energy is ap- 
proximated by 



AS 4 



(17) 



where qo = 1n\J K/k and A an unknown numerical fac- 
tor. The fitting of Eq. Q to the solution of the SCSA 
equations in the region 10~ 4 — 1 A^ 1 gives A = 0.3261, 
as shown in Fig. [4] by the dashed green line. In this fit- 
ting the exponent rj has been fixed to its long wavelength 
value, t) = 0.82. This approximation is a good interpo- 
lation function between the long- and short-wavelength 
regions, and it can be used to simplify the calculation of 
physical quantities that involve the renormalized correla- 
tion function. This range of wavevectors (10 -4 — 1 A" 1 ) 
is relevant for discussing the scattering of electrons by 
ripples 



tor, q c ss 0.1 A -1 , that separates the region of validity of 
the harmonic approximation (for q > q c ) where G(q) oc 
<7~ 4 , from the region where fluctuations lead to a consid- 
erable renormalization of the correlation function, and 
where G(q) oc q~ A+r) . This value of q c is close to the one 
given by the Ginzburg criterion. From this wavevector, 
the exponent rj changes from zero (for q > q c ) to 0.82 in 
the long wavelength limit. This limit is impo rtant w hen 
dealing with MEMS applications of graphene! 9 | 3 ° | 31 1 The 
renormalization of the bending rigidity k —¥ Kji(q) ~ 
should be taken into account, e. g., when calculating the 
eigen-frequencies of graphene membranes that would be- 
come ui(q) ex yj n R (q)q 4 oc q 2 ^^/ 2 ~ q ls . 

Our results show the importance of considering the 
renormalization of the bending rigidity. The good agree- 
ment between SCSA and Monte Carlo simulations for 
graphene can be seen as a proof that SCSA is a good 
approximation to account for the effect of corrugation in 
the physical properties of graphene. 



IV. CONCLUSIONS 

In summary, we have studied numerically the self- 
consistent theory of polymerized membranes proposed in 
Ref. [5D1 The critical exponent that we obtain in the long 
wavelength limit, rj rj 0.82, coincides with the analytic 
approximation. In addition, we have calculated the cor- 
relation function G(q) in the whole range of momenta and 
found good agreement with results of Monte Carlo calcu- 
lations. We have also found the characteristic wavevec- 
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